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In the past decades, Random Electrodynamics (also called Stochastic Electrodynamics) has been 
used to study the classical harmonic oscillator immersed in the classical electromagnetic zero-point 
radiation. Random Electrodynamics (RED) predicts an identical probability distribution for the 
harmonic oscillator compared to the quantum mechanical prediction for the ground state. More- 
over, the Heisenberg minimum uncertainty relation is also recovered with RED. To understand the 
dynamics that gives rise to this probability distribution, we perform an RED simulation and follow 
the motion of the oscillator. This simulation provides insight in the relation between the striking 
different double-peak probability distribution of the classical harmonic oscillator and the Gaussian 
probability distribution of the RED harmonic oscillator. 

A main objective for RED research is to establish to what extent the results of quantum mechanics 
can be obtained. The present simulation method can be applied to other physical systems, and it 
may assist in evaluating the validity range of RED. 

PACS numbers: 03.50.-z, 02.60.-x 



I. INTRODUCTION 

According to quantum electrodynamics, the vacuum 
is not a tranquil place. A background electromagnetic 
field called the zero-point field, or vacuum field, is al- 
ways present, independent of any external electromag- 
netic source [1 . This has inspired an interesting modi- 
fication to classical mechanics, called random electrody- 
namics (RED) [2]. As a variation of classical electrody- 
namics, RED adds a classical zero-point electromagnetic 
field with random phases to the original theory. In this 
paper, the name "RED vacuum field" will be used for 
this classical electromagnetic field. The RED vacuum 
field has some properties identical to the QED vacuum 
fielcQsuch as the Lorentz-invariant spectral energy den- 
sity [H [2|. With the aid of this background field, RED 
is able to reproduce a number of results that were origi- 
nally thought to be pure quantum effects [U [2l [4HZ] • In 
particular, Boyer has shown that the moments (x^) of a 
harmonic oscillator immersed in the RED vacuum field, 
called the RED harmonic oscillator, are identical to those 
of the quantum harmonic oscillator in the ground state 
[8]. As a consequence, the Heisenberg minimum uncer- 
tainty relation is satisfied for an RED harmonic oscilla- 
tor, and the probability distributions of the RED har- 
monic oscillator is also the same as that of the ground 
state quantum harmonic oscillator. 

Although classical mechanics and RED are both theo- 
ries that gives trajectories of particles, the position prob- 
ability distributions of the harmonic oscillator in both 



theories are dramatically different. In classical mechan- 
ics, it is most likely to find the particle at the two turning 
points of the trajectory. Consequently, the probability 
distribution has two peaks as shown in Fig. [l] On the 
other hand, the probability distribution for an RED har- 
monic oscillator is a Gaussian distribution. What hap- 
pens to the dynamics of the classical harmonic oscillator 
in the presence of the RED vacuum field such that its 
probability distributions would become Gaussian? Why 
do the widths of these distributions satisfy Heisenberg's 
minimum uncertainty relation? In this paper, we perform 
an RED simulation that answers these questions. 

Most of the results obtained with the RED the- 
ory reproduce quantum mechanical results analytically. 
Whereas Cole has a series of numerical studies on the 
hydrogen atom 0Hl2], numerical studies in this field are 
rare. The major challenge for RED simulation is to incor- 
porate the RED vacuum field which has an unbounded 
spectrum. A representative selection of vacuum field 
modes is thus the key to successful RED simulations. 
In this study, our approach for the selection of modes 
has been extensively tested for convergence, and we doc- 
ument the details of our simulation so it can be used by 
others. 

The organization of this paper is the following. First, 
in Sec. [IT] Boyer's results on the RED harmonic oscil- 
lator are briefly reviewed. Based on these results, the 
probability distribution for the RED harmonic oscilla- 
tor is derived. Second, in Sec. Illll a detailed numerical 
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method for simulating the RED vacuum field and the 
RED harmonic oscillator is presented. Third, in Sec. IV 
the trajectory of the RED harmonic oscillator is solved 
numerically, and the constructed probability distribution 
is compared to the analytical probability distribution. 
Two sampling methods are used in constructing the prob- 
ability distribution from the simulated trajectories. The 
first approach is "sequential sampling", which is suit- 
able for studying the relation between the dynamics of 
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FIG. 1. A comparison between the classical and the RED harmonic oscillators. Top: Without any external force, the classic 
harmonic oscillator that is initially displaced from equilibrium performs a simple harmonic oscillation with constant oscillation 
amplitude. The resulting probability distribution has peaks at the two turning points. Bottom: Without any external force 
except for the RED vacuum field, the RED harmonic oscillator undergoes a motion that results in a Gaussian probability 
distribution. This motion and its relation to its classical counterpart is investigated in detail in this paper. 



the RED harmonic oscillator and its probability distri- 
bution. The second approach is "ensemble sampling", 
which lends itself well for parallel computing and is con- 



venient for statistical interpretations. Lastly, in Sec. VI 
we discuss some potential applications of RED simulation 
in studying other quantum phenomena, including the 
electron doubleslit diffraction. Numerical studies have 
an advantage over the analytical solutions in that they 
can be adopted to a range of physical systems, and we 
hope that the current method of simulation may assist 
in assessing RED's validity range. 



II. THEORY 



A. Brief Review of Boyer's Work 



In his 1975 papers [218], Boyer calculated the statis- 
tical features of an RED harmonic oscillator, and the 
Heisenberg minimum uncertainty relation is shown to be 
satisfied for such an oscillator. The RED vacuum field 
used in Boyer's work arises from the homogeneous solu- 
tion of Maxwell's equations, which is chosen to be zero 
in classical electrodynamics In unbounded space, the 



RED vacuum field has an integral fornj^ 

2 



A=l ■ 



r?(k,A)E 



(1) 



a(k. A) 



:^(k,A) 



where uj = c|k|, and 6>(k, A) is the random phase uni- 
formly distributed in [0,27r]. The integral is to be taken 
over all k-space. The two unit vectors, s(k, 1) and 
s(k, 2), describe a polarization basis in a plane that is 
perpendicular to the wave vector k. 



s(k. A) •k = 0. 



(2) 



Furthermore, the polarization basis vectors are chosen to 
be mutually orthogonal. 



s(k,l) •s(k,2) = 0. 



(3) 



To investigate the dynamics of the RED harmonic oscil- 
lator, Boyer used the dipole approximation. 



k-r <Cl, 



(4) 



A detailed account of the RED vacuum field in unbounded space 
is given in Appendix I A 11 



3 



to remove the spatial dependence of the RED vacuum 
field, Eq. ([T]), from the equation of motion. Therefore, 
the equation of motion for an RED harmonic oscillator 
used in Boyer's analysis is 

mx = —mujQX + mT'x + eE^^J^{t), (5) 
2e^ 1 
3mc^ 47reo ' 

where m is the mass, e is the charge, uq is the natural 
frequency, and F is the radiation damping parameter. 
The x-component of the RED vacuum field in Eq. ([5| is 



2 

Eiliit) =Y.ld'k s(-)(k, A)^ (a(k, A)e- 

A=l 

+ a*(k,A)e^"^), 
and the steady-state solution is obtained as 



(6) 



7?(k, A) f a(k, A) ^.^^^ 



2 VC(k,A) 



a*(k, A) 



(7) 



where C(k, A) = (— cj^ +^o) — ^Fcj^. Additionally, using 
the condition of sharp resonance. 



(8) 



Boyer further calculated the standard deviation of posi- 
tion and momentum from Eq. ([7| by averaging over the 
random phase |8 , 



2mujQ ' 



(9) 



(10) 



The above result satisfies the Heisenberg minimum un- 
certainty relation. 



h 



(11) 



From an energy argument, Boyer showed that this un- 
certainty relation can be derived from a delicate balance 
between the gain of energy from RED vacuum field driv- 
ing and the loss of energy through radiation damping [2] . 



B. Position Probability Distribution 

Given the knowledge of the moments {x^)q^ the Fourier 
coefficients F^{k) of the position probability distributiorj^ 



^ An RED probability distribution can be constructed from either 
data collected at a sequence of time, or data drawn from an RED 



Pq{x) can be determined through the relation 

/ + 00 
-OO 

= E^/ -"w^- (12) 



n=0 

OO 



n=0 



Using Eq. ^ and the relation from Boyer's paper [8] 

/g±i(^(k,A)+^(k',A'))\ = 0, 



(13) 



the moments {x^)q can be evaluated. 



h 



(2m)! 
m!2™ V 2mwo 



(14) 
(15) 



where m is a natural number. Consequently, only even- 
power terms are contributing in Eq. (12), and the Fourier 
coefficients F^{k) are determined. 



i-ik) 



2m 



m=0 



(2m)! 



(16) 



Therefore, although not explicitly given, it is implied by 
Boyer's work [8 that the position probability distribution 
of the RED harmonic oscillator is 



1 + 



e ^ 



(17) 



which is identical to the position probability distribution 
of the quantum harmonic oscillator in the ground stat^ 



ensemble. In the RED ensemble, each member is characterized 
by the random phase of the RED vacuum field. In this pa- 
per, the probability distribution constructed from the ensemble 
sampling method is denoted as Pq{x). 
^ This result is consistent with the phase space probability distri- 
bution given in Marshall's work [13| . 
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III. NUMERICAL SIMULATION 

A. RED Vacuum Field in Bounded Space 

A field confined in a space of volume V with zero value 
boundary condition has a discrete spectrum, and a sum- 
mation over the wave vectors k is required. On the other 
hand, the field in unbounded space is not subject to any 
boundary condition, so every wave vector k is allowed. 
Thus, the field in unbounded space takes an integral form 
[2 , while the field in bounded space takes a summation 
form [Ills]. In a simulation, it is convenient to write the 
RED vacuum field in its summation form. 



k,A * ^ 



= E 

k,A 



cos(k-r-wt + 6i^ 



(18) 



A^k = k'^ sinOAkAOAcj). To satisfy both conditions, we 
use a set of specifically chosen numbers {f^ijni^ijm^ijn) 
to sample the wave vectors k. Namely, for z = 1 . . . A/"^, 
j = 1 . . . , and n = 1 . . . N^p , 

i^ijn = (^0 - A/2)^/3c^ + - 
^i^n = -1 + (j - I) Ad (22) 
= Rf^ ^{n-l)A^, 



'ijn 



where R^^^ is a random number uniformly distributed in 
[0,27r], and the stepsizes are constant, 



Ak ■ 

Af : 



[{ojQ + A/2)V3c3 - {ojQ - A/2)V3c3 



Wo' 



1' 



(23) 
(24) 
(25) 



Such a set of numbers {i^ijn^^ijn^^ijn) is then used for 
assigning the spherical coordinates to each sampled wave 
vector k. 



where uj = c|k|, 0^^^ is the random phase uniformly dis- 
tributed in [0,27r], and V is the volume of the bounded 
space. A derivation of the summation form of the RED 
vacuum field in bounded space is given in Appendix ] A 2| 
Since the range of the allowed wave vectors k is over 
all k-space and the RED vacuum field diverges for in- 
finitely large |k|, we choose to sample only the wave 
vectors k whose frequencies is within the finite range 
[ooq — A/2,cJo + A/2]. Such sampling is valid as long 
as the chosen frequency range A completely covers the 
characteristic resonance width Tujq of the harmonic os- 
cillator, 



Tuj^ < A. 



(19) 



On the other hand, the distribution of the allowed wave 
vectors k depends on the specific shape of the bounded 
space. In a cubic space of volume the allowed wave 
vectors k are uniformly distributed at cubic grid points, 
and the corresponding RED vacuum field is 



E E 

A=l (k^,ky,k^) 



7^ 



cos{k-r-ujt+0.,)e.,. (20) 



The sampling density is uniform and has a simple relation 
with the space volume F, 



V 



{2n) 



3 • 



(21) 



Nevertheless, such uniform cubic sampling is not con- 
venient for describing a frequency spectrum. In order 
to sample only the wave vectors k in the resonance re- 
gion, spherical coordinates are used. In addition, for 
the sampling to be uniform, each sampled wave vec- 
tor k must occupy the same size of volume element 



where 




kijn^iri^Oijn) COs{(t)ijn) 



kijn^iri (%n)sin 



ki 



, COS ( 



(26) 







U/3 



%jn ~ COS {^ijn) 



(27) 



Consequently, each sampled wave vector k^^^ will be in 
the resonance region and occupy the same size of volume 
element, 

A^k = k^ sinOAkAOAct) 
Ai^Ai^Aif. 

Under the uniform spherical sampling method (described 
by Eqs. (I22|, ([26]), and (p^ expression for the RED 



vacuum field, Eq. (18), becomes 



2 

A=l {K,^,ip) ^ 



cos(k-r-c^t + ^,,Js,,,. (28) 



where k = k^j^. It is worth noticing that when the to- 
tal number of wave vectors A^k becomes very large, both 
uniform spherical and cubic sampling effectively sample 
all the wave vectors k in k-space. Therefore, in the limit 
of large sampling number Ak ^ co, the two sampling 
methods become equivalent]^ and Eq. (21) can be used 



5 The relation pk = = A^k/^k = ^/(27r)3 implies that 
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for both sampling methods to calculate the volume factor 
V in Eqs. (I20l) and (l28|), 



where 



47r/wo + A/2\^ 47r/wo-A/2\^ 



■ 



c J 



(29) 



(30) 



In a computer simulation, the summation indices in 
Eq. (28) can be rewritten as 



2 

A=l i=l j=l n=l 



cos(k-r-a;t + 6'^ Js^^, 
(31) 

where the multiple sums indicate a numerical nested 
loop, and the wave vector k = k^j^ is chosen according to 
the uniform spherical sampling method. In practice, 
and need to be sufficiently large so that the wave vec- 
tor k at a fixed frequency may be sampled isotropically. 
Also, a large is required for representative samplings 
in frequency. As a result, A^k = N^^N^N^ has to be very 
large to achieve convergence when using uniform spheri- 
cal sampling. Therefore, sampling k for many angles at 
each frequency makes RED simulation computationally 
intensive. To improve the efficiency of our simulation, 
we sample k for only one random angle {Oi^(j)i) at each 
frequency (i.e. A^,^ = A"^, A"^ = 1, and A'^^ = !)• Namely, 
A^k N^, and for z = 1 . . . A^^^, 



k, 



where 



ki — 



'ki sin Oi cos ( 
ki sin Oi sin ( 
ki cos Oi 



(3Ki)i/3 
cos-i(t?i) 



(32) 



(33) 



and 



(cjo - A/2f/3c^ + (i - 1)Ak 



R. 



(1) 

i 

(2) 



(34) 



the limit of large sampling number (i.e. — ^ co) is equivalent 
to the limit of unbounded space (i.e. V — oo). At this limit, 
the volume element becomes differential (denoted as dP'k) and is 
free from any specific shape associated with the space boundary. 
Therefore, all sampling methods for the allowed wave vectors 
k become equivalent, and the summation form approaches the 
integral form. This is consistent with the fact that no volume 
factor V is involved in the integral form of the RED vacuum field 
in unbounded space, as shown in Eq. Q. 



The stepsize A/^ is specified in Eq. (23), d = i?*^-^^ is 
a random number uniformly distributed in [—1,1], and 
(f = i?^^^ is another random number uniformly dis- 
tributed in [0, 27r]. As the number of sampled frequencies 
A^^ becomes sufficiently large, the random angles (6>, 0) 
will be uniformly distributed on a unit sphere, and the 
angular distribution of the wave vectors k is isotropic. 

In the limit — <C 1, the above sampling method (de- 



scribed by Eqs. (32), (33), and (34)) and the uniform 
spherical sampling method both approach a uniform sam- 

(jOq 

pling on a spherical surface at the radius Vk = — . In this 



limit, the addition of the condition in Eq. (19) leads to 



the possible choices of the frequency range A, 



Tuo < ^ < 1. 



(35) 



Within this range (Eq. 
vacuum field in Eq. (IS) 



(35|), the expression for the RED 
becomes 



A=l 



cos(k • r - cjt + ^, Jsj^ 



(36) 



where k = k^. For large A^k = 
is calculated using Eq. (29). 



the volume factor V 



Finally, for a complete specification of the simulated 
RED vacuum field, Eq. (36), the polarizations ^ need 
to be chosen. From Eq. (29), we notice that large A'k 
gives large V. Since for large V the RED vacuum field 
is not affected by the space boundary, there is no pref- 
erential polarization direction. Therefore, the polariza- 
tions are isotropically distributed. The construction for 
isotropically distributed polarizations is discussed in de- 
tail in Appendix |B] Here we give the result that satisfies 
the property of isotropy and the properties of polariza- 
tion (described by Eqs. ^ and ([3|), 



Ay) I _ 



f cos Oi cos (j)i cos Xi — sin (^i sin Xi^ 
cos Oi sin (pi cos Xi + cos sin Xi 
— sin Oi cos Xi 




- COS Oi cos (j)i sin Xi — sin (j)i cos Xi\ 

- cos Oi sin (j)i sin Xi + cos (j)i cos Xi , 

sin Oi sin Xi ) 

(37) 



where x is a random number uniformly distributed in 
[0, 2 7r]. With the wave vectors k (described by Eqs. ( [32| ), 
(33), and (34)) and the polarizations e^^ (described by 
Eq. ([37[)), tne endpoints of the sampled RED vacuum 
field vector are plotted on a unit sphere as shown in 
Fig. [2| which illustrates the isotropy of the distribution. 

In summary, the RED vacuum field, Eq. (36), can be 
specified by a set of four numbers (/>:^, i^^, xT), chosen 
according to Eqs. (l32|, ([33]), (|34]), and Eq. The 
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Eq. (18) and the magnetic part B^^c, is 



0.5. 



-0.5. 



-1> 
1 




-1 -1 



FIG. 2. The isotropic distribution of the polarization field 
vectors (number of sampled frequencies number Nco = 
3k). The endpoints of the polarization field vectors, s^^ , are 
plotted for a random sampling of modes (k, 1) according to 
the methods described in Eqs. (32), (33), and (34). 



only assumption used in determining these numbers is 
Eq. (35), which is equivalent to the sharp resonance con- 
dition, Eq. ([8|, used in Boyer's analysis. 



B. Simplified Equation of Motion 

An RED harmonic oscillator is a charged classical har- 
monic oscillator immersed in the RED vacuum field. The 
complete equation of motion is 

mx = -mujlx^mTx E[^^^{x,t) + (v x ^^ac{x,t))^''^ 

(38) 

where m is the mass, e is the charge, and cjq is the 
natural frequency. The radiation damping parameter 
2e2 1 

is r = . Only the motion in the x-direction 

6mc^ 47reo 

is considered in this paper, so the position of the par- 
ticle is r=(x,0,0), and the velocity of the particle 
is v = (x,0,0). Also, the x-component of the vector 

V X B^ac is denoted as (v x B^^c)^^^- 

Three terms are involved in the equation of motion 
— the spring force, radiation damping, and the RED 
vacuum field. 



-F^ spring — UIUJqX^ 
-^damp — X , 

p _ p \^ix) I fv X B )^'^^ 



(39) 
(40) 

(41) 



In unbounded space, the RED vacuum field is given by 
Eq. ([T]). In a bounded space of volume V, the RED 
vacuum field, including both the electric part ^vac as in 



k,A 



^vac{r,t) = ^ J^cos(k-r-cjt + 6'^ Jsj^ 



eoV 



k,A 



B„ac(r,i) = J^cos(k-r-wi + ^^ J(k x e,^ J, 



(42) 



where u = c|k|, and 0^^^ is the random phase uniformly 
distributed in [0,27r]. Throughout this section the RED 
vacuum field as described by Eq. (42) will be used. 



Before a simulation of the RED harmonic oscillator is 
performed, we may increase the efficiency of our com- 
putation by making some approximations to the equa- 
tion of motion, Eq. (38). Namely, the terms Fdamp-, Fvac 
can be approximatedby using the dipole approximation 
k • r <C 1 (Eq. C])), and the sharp resonance condition 
TuJo < 1 (Eq. 

Denoting the order of magnitude for length, time, and 
electric field as xq, to, and Eq respectively, the order 
of magnitude of the two terms of Fyac in Eq. (41) are 
evaluatecrl 



[eE, 



[e(v X B,,,)^-^] 



eEo, 
(^)l 



Cto 



eEo. 



(43) 



The sharp resonance condition Eq. ([8| makes the har- 
monic oscillator responsive to only frequencies that are 
close to ujo, thus the time scale of the particle's motion 

1 Xq 

is to — — . The dimensionless parameter — in Eq. (43) 
UJO cto 
Xq 

can thus be evaluated, — koxo^ to yield 
cto 



[eEill]=eEo, 



[e(v X Byac)^''^] = {koXo)eEo. 



(44) 



Furthermore, the dipole approximation from Eq. (|4| im- 
plies 



koxo < 1, 



(45) 



which can be used to approximate Fyac with its domi- 
nating term, 

Fyac ^ eEo. (46) 
The dipole approximation also removes the spatial de- 



The order of magnitude for dimension A is denoted as [A]. 



7 



pendence in the RED vacuum field, 



E 



k,A 



1 



2! 



(k ■ r)^ 

k-r-^^ + ...)sin(c^t-^,,J 



k,A ^ ^ 



k,A 



(47) 



Through the approximations, Eqs. (46) and (47), the 
equation of motion Eq. (38) becomes 



mx = —mujQX + mT'x + eEl^^{t). 



(48) 



The radiation damping mT'x in Eq. ( |48| ) is weh known to 
give runaway solutions. Boyer's approach is to find the 
solution through Fourier transformation [8 . Following 
this approach, the steady-state solution for Eq. (48) may 
be obtained. 



x{t) 




m 

where C^ ^ = (-cj^+cjg) - iTcu^. Although Eq. (49) 
could be evaluated numerically, this would not allow lor 
straightforward generalization to physical systems where 
the equation of motion is modified, which is one of the 
purposes of this work. Hence, we follow the approach in 
[MHTB] to avoid runaway solutions. According to classi- 
cal electrodynamics, the equation of motion for an elec- 
tron with radiation damping is 



mx = Fext + m^x, 
which can be approximated with 

m^x — Fext 1 



(50) 



(51) 



using the assumption m^'x <C Fe^t- Such an assump- 
tion is necessary for a point-particle description of the 
electron to be valid [TT, T6^. With the reduced equa- 
tion (Eq. (51)), the radiation damping m^'x may be es- 
timated as 

m-fx :^ -fPext, (52) 
and then iterated back to the original equation (Eq. ([50|), 

mx Fext ^iFext' (53) 

This approximate equation of motion is free of runaway 
solutions. Applying Eq. (53) to the RED harmonic oscil- 
lator, Eq. (I48|), we obtain 



mx 



-mujnX - 



■mTco^x + eEill + eTEill (54) 



The order of magnitude of each force term is 



F 

spring 
vac 



ImujnXl 



(^)i 



: mujQXo^ 
eEo, 



F, 



damp 



[mVuoQx] = (rco'o)mcjo^o 



(55) 



To evaluate the magnitude of the terms in Fdampi a 
random walk model may be used. For a fixed time 
t = to^ a particular value of Ev^c{to) and x{to) are found 
in Eqs. (|47| and (|49|, and the order of magnitude for 

^iac(^o) and x{to) are equal to and xq. The mathe- 
matical form of the complex £'iac(to) and x{to) is analo- 
gous to a two-dimensional random walk on the complex 
plane with random variable Ojk^Aj^ where {k. A} denotes 
a set of modes (k. A). Thus, if averaged over {k. A}, the 
order of magnitude for £^iac(^o) and x{to) may be esti- 
mated by the root-mean-squared distance of ^iac(^o) and 
x{to). In a two-dimensional random walk model [17 , the 
root-mean-squared distance Drms is given by 

Drms = VN~s ■ As, (56) 

where Ng is the number of steps taken, and As is a typ- 



ical stepsize; for Drms-, 
e 



As 



\\ ~T~r^ and for 



As 



1 



2 \mVuol \l cqV I' 
tude, £^0 and xq, may be estimated a^ 



Hence, the order of magni- 



Xq D, 



ix) 




The order of magnitude for F^amp 
ingly. 



(57) 



is evaluated accord- 



[mFwoi] ~ e 



2 V eoV 



eoV 



(58) 



Again, using the sharp resonance condition Tuq <C 1 
(Eq. ([8|), Fdamp is approximated with its dominating 
term. 



Fdamp — 6 



I hup 



(59) 



Using V = {27t)^Nu;/Vi^ and - 47rcjg (Fcjg) /c^ (Eqs. ( |29| > 
a nd (|30|)), the value of xq can be estimated as xq ~ 
y^3/7r h/2mujo , which is consistent with Boyer's calculation for 
the standard deviation of position in Eq. joj). 
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and the equation of motion, Eq. (54), is further simphfied 

tcE] 



mx :^ —muj. 



ix-mTojlx + eEiliit). (60) 



As an additional note, given the estimation of Eq and 
xq in Eq. (57), the three force terms in Eq. (60) have the 



fohowing relation. 



\mujnX 



spring 



(61) 



domi- 



The implication is that the spring force F< 
nates, while the RED vacuum field Fyac and the radia- 
tion damping Fdamp are perturbations. The driving force 
Fyac and the damping force Fdamp will establish a bal- 
ance, keeping the amplitude of the motion of the RED 
harmonic oscillator in the vicinity of xq. 

Finally, after the main equations for simulation are 
completely defined (Eqs. (47) and (60)), the range of the 
integration time Tint (i-e. how long the simulation runs) 
needs to be chosen. Upon careful inspection, two im- 
portant time scales can be identified from the analytical 
solution of Eq. (60), 



m 



k,A 



^1 ( Si^^-iujt 



■k,A 



k,A ' 



satisfy Tint < ^rep to avoid repetitive solutions. A de- 
tailed discussion about Trep can be found in Appendix [Cj 
Here we give the choice of r^^^. 



27r 



(67) 



where Auj is the gap between adjacent frequencies. The 
frequency gap Au can be estimated using Eqs. (23), (33), 
and dMl), 



3 K.0 ' 



(68) 



where = - f — ^ . 

3 V c / 

To summarize, Eq. ( [6Q| ) is the equation of motion that 
will be used for simulating the RED harmonic oscilla- 
tor. The RED vacuum field in Eq. (60) is simulated by 



Eq. (47), and the specifications of its modes (k. A), polar- 
izationss^, and other relevant variables can be found 
in SeclIII A] The simplification from Eq. (ISSl) to Eq. (leol) 



relies on only two conditions — namely, the dipole ap- 
proximation Eq. (|4| and the sharp resonance condition 
Eq. ([8|, which are also the only two conditions used in 
Boyer's analysis [8 . The parameters e,m, and uq should 
be chosen to satisfy these two conditions. Lastly, the 
integration time Tint of the simulation is to be chosen 
within the range nrans < nm < ^rep, where Ttran and 



(62) Trep are given in Eqs. (66) and (67) respectively. 



where A is a coefficient determined by the initial condi- 
tions, and 



IV. RESULTS 




(63) 

(64) 
(65) 



The first term in Eq. (62) is the transient part of the 



solution. The characteristic time for the transient motion 
is 



Ttr 



(66) 



The position should be evaluated for Ttrans <^ '^int^ if 
we are only interested in the steady-state motion. The 
second term in Eq. (62) is the steady-state part of the 
solution. Because the steady-state solution is a finite 
discrete sum of periodic functions, it would have a non- 
physical repetition time Trep- The choice of Tint should 



Both Boyer's equation of motion, Eq. ([s]), and the equation of 
motion used in our work, Eq. ( |60| ), give the same result on the 
standard deviation (Eq. Jol) and the moments (Eq. ITsl). 



In Sec. |IIB[ it was shown that the RED harmonic os- 
cillator probability distribution is Gaussian, contrasting 
to the double-peak probability distribution of the clas- 
sical harmonic oscillator. We set out to investigate how 
the dynamics of the RED harmonic oscillator can give 
rise to such a Gaussian distribution. In this section, the 
result of the simulation for an RED harmonic oscillator 
is presented, and the relation between the trajectory and 
the probability distribution is discussed. 

In constructing the probability distribution from the 
particle trajectory, two sampling methods are used, 
namely sequential sampling and ensemble sampling. In a 
sequential sampling, the position or velocity is recorded 
in a time sequence from one trajectory, while in an en- 
semble sampling, the same is recored only at the end 
of the simulation from an ensemble of trajectories. The 
recorded positions or velocities are collected in histogram 
and then converted to a probability distribution, which 
is compared to the analytical result, Eq. ( p!7| ). Whereas 
the sequential sampling illustrates the relationship be- 
tween the buildup of the probability distribution and the 
dynamics of the trajectory, the ensemble sampling is con- 
venient for statistical interpretation. In addition, the en- 
semble sampling is suitable for parallel computing, which 
improves the computation efficiency. 
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FIG. 3. A comparison between the field evolution and the 
particle trajectory. Top: RED vacuum field (red solid line) 
is compared to the trajectory of the RED harmonic oscilla- 
tor (black solid line). Bottom: A magnified section of the 
trajectory shows that there is not a fixed phase or amplitude 
relation between the particle trajectory (black solid line) and 
the instantaneous driving field (red solid line). The modu- 
lation time for the field (red solid line) is also shown to be 
longer than that for the motion of the harmonic oscillator 
(black solid line). 




Momentunn(kg-nn/s) xlO"^^ Time(s) ^ ^q-i2 



FIG. 4. The RED probability distribution constructed from 
a sequential sampling (number of sampled frequencies Nco = 
20k). Left: The position and momentum probability distri- 
butions for the RED harmonic oscillator (black dot) and the 
ground state quantum harmonic oscillator (red & blue solid 
line) are compared. Right: A illustration of the sequential 
sampling shows how positions are recorded at a regular time 
sequence (red cross). Note that at small time scale the oscil- 
lation amplitude is constant, but at large time scale it modu- 
lates. 



A. Sequential Sampling and the Relation to its 
Probability Distribution 



By solving Eq. ( |6Q[ ) numerically, the steady-state tra- 
jectory for the RED harmonic oscillator is obtained as 
shown in Fig. [3] For a comparison, the RED vacuum 
field is also included. The differential equation (Eq. ([60|) 
is solved using the adaptive 5th order Cash-Karp Runge- 
Kutta method [TF, and the integration stepsize is set 

as one twentieth of the natural period, — ( — |. The 

charge e, mass m, and natural frequency ujq are chosen 
as (7e, 10~^me, and 10^^ (rad/s) respectively, where Qe 
is the electron charge and mg is the electron mass. The 
frequency range A of the RED vacuum field is chosen to 
be much broader than the resonance width, 220 X Tuj^ 
The choice of 10~^me is made to bring the modulation 
time and the natural period of the harmonic oscillator 
closer to each other. In other words, the differential equa- 
tion (Eq. ([60])) covers time scales at two extremes, and 
the choice of mass 10~^me brings these scales closer to 
each other while keeping the integration time manageable 
without losing the physical characteristics of the prob- 
lem. Some features of the trajectory, shown in Fig. |3j 
are worth mentioning. First, there is not a fixed phase 
or amplitude relation between the particle trajectory and 
the instantaneous driving field. Second, the rate of am- 
plitude modulation in the particle trajectory is slower 
than that in the driving field. While these results may 



seem a little surprising at first glance, the steady-state 
solution in Eq. (62) provides some insights into both of 
these behaviors. 

To understand the features of the trajectory, we study 
the temporal response of the harmonic oscillator as re- 
vealed in the Green formulation of the steady-state solu- 
tion [H], 



x{t) 



mujR 



J — oo 



-r^o'(*-t')/2 gij^ ^ajRit-t'))dt'. 

(69) 

The solution indicates that the driving field E^^^{t^) at 
any given time f has an effect on the particle for a finite 

period of time, ^r^- other words, the particle trajec- 
tory x{t) at a particular time t is affected by the driving 
field E^^^{t') from all previous moments t' <t. This ex- 
plains why the particle trajectory does not seem to have 
a fixed phase and amplitude relation with the instanta- 
neous driving field. Another implication of Eq. (69) is 

that it always takes a finite period of time, ip^-^, for the 

particle to dissipate the energy gained from the driving 
field. Thus, even if the field already changes its am- 
plitude, it would take a while before the particle can re- 
spond. This is why the amplitude modulation in the par- 
ticle trajectory is slower compared to that in the driving 
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FIG. 5. Contributions of different oscillation amplitudes in 
the final probability distribution. Top: A steady-state trajec- 
tory of the RED harmonic oscillator (black) and its several 
sections (red) are shown. A section is limited to the dura- 
tion of a modulation time. The oscillation amplitude changes 
significantly outside a modulation time, so different sections 
have different oscillation amplitude. Bottom: The probability 
distributions of each trajectory section. Since the oscillation 
amplitude is approximately constant in each section, the cor- 
responding probability distribution (red bar) is close to the 
classical double-peak distribution. These probability distri- 
butions contribute to different areas of the final probability 
distribution (black dash line) that is constructed from the full 
steady- state trajectory. 

fiel(fl 

Using a sequential sampling method, a probability dis- 
tribution is constructed from the simulated trajectory 
of a RED harmonic oscillator as shown in Fig. |4] The 
RED probability distribution turns out to be a Gaus- 
sian distribution, which is the same as the ground state 
probability distribution of a quantum harmonic oscilla- 
tor. To understand the relation between the trajectory 
and the Gaussian probability distribution, we investigate 
the dynamics of the RED harmonic oscillator at two time 
scales. First, at small time scale, the particle oscillates 
as a classical harmonic oscillator. The oscillation ampli- 
tude is constant, and the period is T = 27t/ujo. Such 
an oscillation makes a classical double-peak probability 
distribution. Second, at large time scale, the oscillation 
amplitude modulates. As a result, different parts of the 
trajectory have double-peak probability distributions as- 
sociated with different oscillation amplitudes, which add 
to make the final probability distribution a Gaussian dis- 
tribution as shown in Fig. |5] 



^ However, in case of slow field modulation (i.e. field bandwidth 
shorter than resonance width of the harmonic oscillator) the 
modulation time of the field and the particle are the same. 




FIG. 6. A schematic illustration of the oscillation amplitude 
as a sum of different frequency components in the complex 
plane. At a particular time t = to, an oscillation amplitude 
A (to) (blue solid arrow) is formed by a group of frequency 
components (blue dash arrow), which rotate in the complex 
plane at different rate. After a time At ^ Tcoh, the angles 
of the frequency components (red dash arrow) change only 
a little. Therefore, the oscillation amplitude A{to + At) (red 
solid arrow) does not change much within the coherence time. 



To further understand how amplitude modulation 
turns the final probability distribution into a Gaussian 
distribution, we again inspect the steady-state solution 
in Eq. (|62|), 



E 

k,A 



huj 1 

^2 \C, 



C.C. 



Ax) 



(70) 



The solution can be rewritten in terms of the complex 
position x(t). 



x{t) = Re (x(t)). 



(71) 



Since the frequencies are sampled symmetrically around 
cjo, the amplitude modulation is made evident as the 
complex position x{t) is factorized into a modulation 
term A{t) and an oscillation term 6"*^^°^, 



x(t) = A(t)e-^^°S 
A{t) 



k,A 



-i{oj—ojo)t 



(72) 



The amplitude A{t) is a sum of many complex compo- 
nents, which rotate in the complex plane at different 
rates, u — uq. At any given time, the amplitude A{t) 
is determined by the configuration of these components, 
as shown in Fig. |6] As time elapses, the configuration 



11 




FIG. 7. Sampled oscillation amplitude and the reconstructed position probability distribution. Top: The oscillation amplitudes 
(red dot) are sampled from the full steady- state particle trajectory (black) with a sampling time- step larger than one coherence 
time Tcoh- Middle: The sampling time-step for the oscillation amplitudes (red dot) is approximately Srcoh- Bottom: The 
probability distributions of the position (black dot) and the oscillation amplitude (red bar) are constructed from the simulated 
data. Notice that the most frequently occurring oscillation amplitude corresponds to the half-maximum width of the position 
distribution (black dot), which is offset by As = 1.5 x 10^. The reconstructed position probability distribution (blue solid line) is 

given by P{x) = Pa{x) = / PA{x)f{A)dA, where Pa{x) is the classical double-peak probability distribution associated with 

A 

an oscillation amplitude A, and f{A) is the probability distribution of the oscillation amplitude. This procedure demonstrates 
that the Gaussian RED probability distribution can be thought of as a sum of classical double-peak distributions with a specific 
amplitude distribution, answering the question illustrated in Fig. ^ 



of the components evolves and the amplitude A{t) may 
change dramatically. However, when the elapsed time 
At is much smaller than the coherence time Tco/i, which 
is given by the shortest rotating period, the change of 
amplitude A{t) is negligible. 



A{t + At) 



^coh 



A(t) for 
27r 
_\uj - CJol 



At < Tcoh, 



(73) 
(74) 



Therefore, since an oscillation amplitude A{t) is approx- 
imately constant within one coherence time Tcq/J^ the 



The coherence time TcoK defined in this work has the same mean- 
ing as the temporal width of the first-order correlation function. 
We can take the autocorrelation of the simulated trajectory and 
show that its temporal width is the same as the coherence time 
computed in this work. 



amplitude modulation as previously discussed can be 
characterized by the coherence time Tcoh- In other words, 
the coherence time Tcoh is the modulation time. 

Now, if we take a representative sampling of the oscil- 
lation amplitudes (i.e. each sampled amplitude is several 
coherence times Tcoh apart), we obtain an amplitude dis- 
tribution as shown in Fig. [7| From the amplitude dis- 
tribution, it is clear that the occurrence of very large 
or very small amplitudes is rare. This is because it re- 
quires all the components of the amplitude A{t) to be 
either completely aligned or completely misaligned to ob- 
tain the extreme values of A(t). For most of the time, 
the configuration is in partial alignment, and the ampli- 
tude A{t) is of the average size as calculated previously 
(i.e. xq). Using the amplitude distribution in Fig. [7| a 
Gaussian distribution can be reconstructed by adding up 
all the classical double-peak distributions associated with 
these sampled oscillation amplitudes. This reconstructed 
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FIG. 8. The RED probability distribution constructed from 
an ensemble sampling (number of particles Np — 200k, num- 
ber of sampled frequencies Nco = 2k) in contrast to the se- 
quential sampling (Fig. [4]). Left: The position and momen- 
tum probability distributions are shown for the RED har- 
monic oscillator (black dot) and the ground state quantum 
harmonic oscillator (red & blue solid line). As the ensem- 
ble sampling corresponds to the procedure of random phase 
averagmg (Eqs. ^ and ([lO}), the widths of the two RED dis- 
tributions satisfy Heisenberg's minimum uncertainty relation, 
cTxCTp = h/2. Right: A illustration of the ensemble sampling 
shows how positions (red cross) are recorded from an ensem- 
ble of trajectories (black). 



Gaussian probability distribution resembles the proba- 
bility distribution obtained directly from the simulated 
data as shown in Fig. ^ It is worth pointing out that 
it is impossible to reconstruct such a Gaussian distri- 
bution with any amplitude distribution other than the 
one given by the particle trajectory (Fig. [t]). This result 
shows that the Gaussian probability distribution of the 
RED harmonic oscillator is a consequence of its intrinsic 
amplitude modulation. 



B. Ensemble Sampling and Heisenberg 
Uncertainty Relation 

In many RED analyses [H [HI [20], the procedure 
of random phase averaging is an essential tool in obtain- 
ing the statistical properties of a physical system. To 
compare the RED numerical results with the analyses, it 
is advantageous to use ensemble sampling. In ensemble 
sampling, particles in each trajectory of the ensemble are 
prepared with identical initial conditions, but the RED 
vacuum fields differ in the initial random phases 6^^^. 
Here, different random initial phases 0^ ^ represent the 
physical realization of random phase averaging. At the 
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FIG. 9. Radiation damping and Heisenberg's minimum un- 
certainty relation. Under the balance between vacuum field 
driving and the radiation damping, the trajectory of the RED 
harmonic oscillator (black solid line) satisfies the Heisenberg 
minimum uncertainty relation. When the radiation damping 
is turned off in the simulation (red & blue solid line) , the mini- 
mum uncertainty relation no longer holds, although the range 
of the particle's motion is still constrained by the harmonic 
potential. 



end of the simulation, physical quantities such as position 
and momentum are recorded from an ensemble of trajec- 
tories as shown in Fig. [S] In the case of RED harmonic 
oscillator, the recorded positions and momentum are con- 
verted into probability distributions, and the Heisenberg 
minimum uncertainty relation is found to be satisfied. 
Boyer suggested that the minimum uncertainty relation 
is established by a delicate balance between the gain of 
energy from RED vacuum field driving and the loss of en- 
ergy through radiation damping [2 . As we turn off the 
radiation damping in the simulation (Fig. [9|, the min- 
imum uncertainty relation no long holds, although the 
range of the particle's motion is still constrained by the 
harmonic potential. This result supports the balancing 
mechanism proposed by Boyer. 

Unlike sequential sampling, ensemble sampling has the 
advantage that the recorded data are fully uncorrelated. 
As a result, the integration time does not need to be 
much longer than one coherence time Tcoh after the parti- 
cle reaches its steady-state. However, this does not mean 
that the computation time will be shorter. In fact, be- 
cause only one data point is recorded from each trajec- 
tory, a simulation with ensemble sampling usually takes 
much longer than with sequential sampling. For exam- 
ple, a typical RED simulation with sequential sampling 
(i.e. number of sampled frequencies = 20k) takes 2.3 
hours to finish, but a typical RED simulation with ensem- 
ble sampling (i.e. number of particles Np = 200/c, number 
of sampled frequencies = 0.5k) takes 61 hours, which 
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FIG. 10. The inverse relation between the computation time 
and the number of processors. The data (black open square, 
blue & red dot) follow an inverse relation y = C/x (black, 
blue & red dash line). The parameter C is determined by 
the product of Np (number of particles) and Nuj (number of 
sampled frequencies). 



FIG. 11. The ensemble-averaged energy of the RED harmonic 
oscillator and its convergence as a function Nuj. The number 
of particles is Np = 200k, and the sampled frequency range 
is A = 220 X Tuoq. The energy converges as the number 
of sampled frequencies N^j increases. At N^j = 0.5k (blue 
cross) , the deviation from the ground state energy of quantum 
harmonic oscillator is 1%. 



is 27 times longer. 

A resolution to this problem is to incorporate parallel 
computing in the ensemble sampling approach. The par- 
allelization of the RED simulation is straightforward for 
ensemble sampling, since each trajectory are independent 
except for the random initial phases ^ . To reduce the 
amount of interprocessor communication and computa- 
tion overhead, each processor is assigned an equal amount 
of work. A typical simulation with ensemble sampling ap- 
proach (i.e. number of particles Np = 200k, number of 
sampled frequencies N^ = 0.5k) takes 61 hours to run 
on a single processor. After parallelization, it takes only 
4.5 hours to run on 80 processors, which is a factor of 14 
speed-upp^ and the computation time scales inversely 
with the number of processors as shown in Fig. [Toj In 
addition, the parallelized code is advantageous for testing 
the numerical convergence of the simulation. In Fig. [TT] 
the convergence of the ensemble- averaged energy as a 
function of sampled frequency number is shown. 



V. SUMMARY & CONCLUSIONS 

The analytical probability distribution of an RED har- 
monic oscillator is obtained in Sec. IIIBI The numerical 
method for the RED simulation is documented in Sec. lIIII 
Agreement is found between the simulation and the an- 
alytical results. Namely, the probability distributions 
constructed from the simulated data follows the same 
Gaussian distributions as the ground state quantum har- 
monic oscillator (Figs. [4j and[8|. As a consequence, the 
Heisenberg minimum uncertainty relation is satisfied by 
the widths of the RED position and momentum distri- 
butions. In the absence of radiation damping, this min- 
imum uncertain relation no longer holds (Fig. [9|, which 
supports Boyer's energy argument that there is a deli- 
cate balance between the vacuum field driving and the 
radiation damping. Overall, the RED harmonic oscil- 
lator is distinct from the classical harmonic oscillator in 
that its oscillation amplitude modulates at the time scale 
of coherence time Tcoh (Fig. 12), and the resulting RED 
probability distribution can be thought of as a sum of 
classical double-peak distributions with a specific ampli- 
tude distribution (Fig. [t]). 



The parallelization of the RED simulation program is devel- 
oped and benchmarked with assistance from University of Ne- 
braksa Holland Computing Center. The program is written in 
Fortran and parallelized using Message Passing Interface (MPI) 
[TSlET]. The compiler used in this work is from the GNU Com- 
piler Collection (GCC) "openmpi-1.4.3/gcc-4.4. 1". However, 
with the commercial compiler "openmpi-1 .4.4/intel-ll . 1" , we 
have found that the computation time can be further reduced by 
a factor of 6. 



VI. DISCUSSION^APPLICATION OF RED 
SIMULATION TO OTHER PHYSICAL SYSTEMS 

In quantum mechanics, the harmonic oscillator has ex- 
cited, coherent, and squeezed states. A natural extension 
of our current work is to search for the RED correspon- 
dence of such states. Currently, we are investigating how 
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FIG. 12. A comparison between the classical and the RED harmonic oscillators. Top: A free classic harmonic oscillator that 
is initially displaced from equilibrium demonstrates a motion of constant amplitude, which results in double-peak probability 
distribution. Bottom: Compared to Fig. [l] it is clear now that the RED harmonic oscillator undergoes an oscillatory motion 
with modulating oscillation amplitude. The oscillation amplitude modulates at the time scale of coherence time TcoH and is 
responsible for the resulting Gaussian probability distribution. 



a Gaussian pulse with different harmonics of uoq will af- 
fect the RED harmonic oscillator. Can the RED har- 
monic oscillator support a discrete absorption spectrum, 
and if so, how does it compare with the prediction from 
quantum mechanics? Such a study is interesting in the 
broader view of Milonni's comment that RED is unable 
to account for the discrete energy levels of interacting 
atoms P, and Boyer's comment that at present the line 
spectra of atoms is still unexplained in RED [22 . 

The current method of RED simulation may also be 
applicable to study other quantum systems that are re- 
lated to the harmonic oscillator, such as a charged parti- 
cle in a uniform magnetic field and the anharmonic oscil- 
lator [4l[20]. For the first example, classically, a particle 
in a uniform magnetic field performs cyclotron motion. 
Such a system can be viewed as a two-dimensional os- 
cillator, having the natural frequency set by the Larmor 
frequency. On the other hand, a quantum mechanical 
calculation for the same system reveals Landau quanti- 
zation. The quantum orbitals of cyclotron motion are 
discrete and degenerate. Such a system presents a chal- 
lenge to RED. For the second example, a harmonic po- 
tential can be modified to include anharmonic terms of 
various strength. Heisenberg considered such a system 
a critical test in the early development of quantum me- 
chanics [23 , 24 . We think that a study of the anharmonic 
oscillator is thus a natural extension of our current study 
and may serve as a test for RED. 

Lastly, over the last decades there has been a sus- 
tained interest to explain the origin of electron spin and 
the mechanism behind the electron doubleslit diffraction 
with RED ^25rf28j . Several attempts were made to con- 



struct a dynamical model that accounts for electron spin. 
In 1982, de la Pefia calculated the phase averaged me- 
chanical angular momentum of a three-dimensional har- 
monic oscillator. The result deviates from the electron 
spin magnitude by a factor of 2 [2F. One year later, 
Sachidanandam derived the intrinsic spin one-half of a 
free electron in a uniform magnetic field [26]. Whereas 
Sachidanandam's calculation is based on the phase av- 
eraged canonical angular momentum, his result is con- 
sistent with Boyer's earlier work where Landau diamag- 
netism is derived via the phase averaged mechanical an- 
gular momentum of a free electron in a uniform magnetic 
field [4 . Although these results are intriguing, the most 
important aspect of spin, the spin quantization, has not 
been shown. If passed through a Stern-Gerlach magnet, 
will the electrons in the RED description split into two 
groups of trajectoriej^ At this point, the dynamics 
becomes delicate and rather complex. To further inves- 
tigate such a model of spin, an RED simulation may be 
helpful. 

On the other hand, over the years claims have been 
made that RED may predict doubleslit electron diffrac- 
tion [2l[27l[28]. In order to explain the experimentally ob- 
served electron doubleslit diffractior|^ [ 351(37] , different 



^ Electron Stern-Gerlach effect is an interesting but controversial 
topic in its own right. Whereas Bohr and Pauli asserted that an 
electron beam can not be separated by spin based on the concept 
of classical trajectories 29 , Batelaan and Dehmelt argue that 
one can do so with certain Stern-Gerlach-like devices 30-33 . 

^ See 34j for a movie of the single electron buildup of a doubleslit 
diffraction pattern. 
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FIG. 13. Schematic illustration of RED doubleslit diffraction. 
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Appendix A: RED Vacuum Field in Unbounded and 
Bounded Space 



In "unbounded" space, the modes are continuous 
and the field is expressed in terms of an integral. In 
"bounded" space, the modes are discrete and the field is 
expressed in terms of a summation. In both cases, the 
expression fo r the field amplitude needs to be obtained 
(see Al and A 2). The integral expression helps com- 



parison with analytical calculations in previous papers 
[H IIHBI E] , while the summation expression is what we 
use in our numerical work. 



1. Unbounded Space 



mechanisms motivated by RED were proposed [27l [28] , 
while a detailed calculation of the diffracting vacuum field 
in the vicinity of a slit structure has been given [38 . In 
1999, Kracklauer suggested that particles steered by the 
modulating waves of the RED vacuum field should dis- 
play a diffraction pattern when passing through a slit, 
since the RED vacuum field itself is diffracted [27 . In 
recent years, another diffraction mechanism is proposed 
by Cavalleri et al. in relation to a postulated electron 
spin motion ^Hj. Despite of these efforts, Boyer points 
out in a recent review article that at present there is still 
not a concrete RED calculation on the doubleslit diffrac- 
tion ^2 . Boyer suggests that as the correlation function 
of the RED vacuum field near the slits is modified by 
the slit boundary, the motion of the electron near the 
slits should be influenced as well. Can the scattering of 
the RED vacuum field be the central mechanism behind 



the electron doubleslit diffraction (Fig. 13)? As Heisen- 
berg's uncertainty relation is a central feature in all mat- 
ter diffraction phenomena, any proposed mechanism for 
electron doubleslit diffraction must be able to account for 
Heisenberg's uncertainty relation. This is exactly what 
RED does for the harmonic oscillator. Thus, we hope 
that the current simulation method may help providing 
a detailed investigation for any proposed RED diffraction 
mechanism. 



The homogeneous solution of Maxwell's equations in 
unbounded space is equivalent to the solution for a wave 
equation. 



E(r, *) =^ E / d'k e(k, A) (i(k, A)e^('^-- 
+A*(k, A)e-*('^""-"*)) 



B(r, t)=-J2 d'k (k X £(k, A)) (i(k, X)e' 

^ A=l 



(k-r-ujt) 



+A*{k, A)e-*('^--"*)) , 



(Al) 



where A(k, A) is the undetermined field amplitude for the 
mode (k. A) and has the unit of electric field (V/m), k 
is defined as the unit vector of k, and the two vectors, 
s(k, 1) and s(k, 2), describe an orthonormal polarization 
basis in a plane that is perpendicular to the wave vector 
k. 

Without loss of generality, a random phase 
g*6'(k,A) -^^ factored out from the field amplitude 
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A(k, A) = A(k, A)e 

2 



ie(k,A) 



E(r, t)=^Y. f"^'^ ^(^' ^) (^(^' ^)^' 

A=l 



(kT-u;^)gi^(k,A) 



+A*(k,A)e- 



-i(kT-u;^)g-i6>(k,A) 



1 ^ r 

^ A=l 



gi(k.r-u;t)gi0(k,A) 



(A2) 



where |(k, A) = k x s(k, A). The field amphtude A(k, A) 
can be determined through the phase averaged energy 
density, 



2mo 



(A3) 



where eo is the vacuum permittivity, /io is the vac uum 
permeabihty, and is the random phase in Eq. (A2). To 
evaluate the phase averaged energy density (li)^, we first 
calculate |E|%nd |B|^ using Eq. ( [A2| , 

|E|'(r,t) = E(r,t)E*(r,t) 

= ZE fd'kfd'k^ (s(k,A)•s(k^AO)/, 

|B|'(r,t) = B(r,t)B*(r,t) 

= A'^A'^' (|(k,A)•|(k^AO)/, 

(A4) 



where / = /(k, A; k', A'; r, t), 
/(k,A;k',A';r,t) 

= Aik, X)A*{k', A')e'(k-r--t+e"(k,A))g-i(k'.r-c.'t+9(k',A')) 
+ A(k, A)A(k', X')e'^^-''-"t+9{k,X))^i(k'-r-u.'t+9{k',X')) 



Consequently, ^|E| ^y. and ^|B|^y . can be evaluated us- 
ing Eqs. (|A4| and (AJ]), 



El' 



lE A'fcA'fc' (£(k,A)-£(k',A'))(/),- 

A=l 

-^^ /d^Jd^fc' (^(k,A).C(k',A')){/),- 
^.jzfd'k \A{k,X)\\ 



(A8) 

The above calculation leads to a relation between the 
field amplitude A(k, A) and the phase averaged energy 
density {u)q in unbounded space. 



A=l 



(A9) 



Now, if we postulate that the vacuum energy is huj/2 
for each mode (k. A), then in a bounded cubic space of 
volume V the vacuum energy density is 

1 



2 

k,A 

1 ( 

k,A ^ 



Akx 



27r 



Ak, 



A=l ^ ^ k 



(AlO) 



In the limit of unbounded space (i.e. V oo), the volume 
A*(k A)74*(k' y^g-i(kT-u;^+^(k,A))g-i(k'T-u;'^+^(k',A')) element A^k becomes differential (i.e. A^k d?k) and 

A*(k, A)A(k^ y)e-^^^-'^-"^+^>'^))e^^^'-'^-"'^+^>''^')). 

t^^^^l- (An. 



the vacuum energy density becomes 

2 



The random phase average can be calculated with the 
following relation [8], 



A=l 



g±i(^~(k,A)+0>',A'))\ = 0, 



Comparing this result with Eq. (A9), we find 



(A6) 



>±i(e(k,A)-e(k',A')) 



5A',A<53(k' - k). 



|^.ac(k,A)r 



Tiuj 



(27r)3eo' 



(A12) 



Applying Eq. ( |A6[ ) to Eq. ( |A5[ ), we obtain 
(/(k,A;k^y;r,t)),~ 

= A(k, A)A*(k^ y)e^^^-^"^'^e"^^^'-^"^''^(5A',A^^(k' - k) 
+ A*(k, A)A(k^ y)e"^^^-^"^'^e^^^'-^"^''^(^A',A^^^(k' - k). 



Assuming A^ac(k, A) is a positive real number, the vac- 
uum field amplitude in unbounded space is determined. 



^i;ac(k. A) — 



flUO 



(A13) 



(A7) Therefore, the RED vacuum field in unbounded space is 
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found to be 



gi(k.r-u;t)gi^(k,A) 



A=l 

^g-i(k-r-u;^)g-i^(k 



A=l 



i(kT-u;t)g-i6'(k,A) 



(A14) 



which is Eq. ([T]). 



2. Bounded Space 



Comparing Eq. (A17) and Eq. (A18), the vacuum field 



amphtude in bounded space is determined, 

^.ac, = (A19) 

Therefore, the RED vacuum field in bounded space is 



2 V V eoV- 

k,A \ 



k,A 



(A20) 



The solution of homogeneous Maxweh's equations in 
bounded space has the summation form, 



k,A 
k,A 

(A15) 

where ^ = k x ^ , ^ is the undetermined field 
amplitude for the mode (k. A) and has the unit of electric 
field (V/m), k is defined as the unit vector of k, and 
the two vectors, Sj^^ and e^^^, describe an orthonormal 
polarization basis in a plane that is perpendicular to the 
wave vector k. 

Using the relatiorj^ 







(A16) 



k',k, 



we can follow the same argument in A 1 and obtain the 
phase averaged energy density in bounded space 



IE 1^: 



2 

k,A I ' 



(A17) 



k,A 



where A^^ ^ = A^^e . Again, if we postulate that the 
vacuum energy is fict;/2 for each mode (k. A), then in a 
bounded space of volume V the vacuum energy density 
is 



1 



(A18) 



k,A 



If the two modes are not identical (i.e. \<i' ^ \<i or X' ^ X) ^ then 



e*^^'^ and e*^^''^' are independent, which leads to the factoriza- 
tion {e^(^k,A±V,A')^_ _ (e^ek,A)-(e±^^k',A')- = q. 



Appendix B: Isotropic Polarization Vectors 

A wave vector chosen along the z-axis, 

/O 


has an orthonormal polarization basis in the xy-plane, 



k = U J = I 



(Bl) 





f COS X] 




(-sin x\ 








cosx 








^ / 



(B2) 



where the random angle x is uniformly distributed in 
[0,27r]. To obtain the wave vector k in Eq. (32), k can 
be first rotated counterclockwise about the y-axis by an 
angle 6>, then counterclockwise about the z-axis by an 
angle (j). The corresponding rotation matrix is described 
by 



R = R):^R 



'cos(j) —sin(j) 0^ 
sm (j) cos (j) 

1; 

' cos cos (j) — sm ( 
cos sin (j) cos ( 
- sin 6> 



and k is obtained accordingly, 

sin cos ( 
k = Rk = I ksinO sin c 
k cos 



cos sin i 

1 
— sin 6> cos ( 



(B3) 




(B4) 



In the same manner, we can rotate s, , and s, ^ with 

" k, 1 k,2 

the rotation matrix i?, and obtain an isotropically dis- 
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FIG. 14. Repetition time of a beat wave. A beat wave (black 
solid line) is made of two frequency components (red & blue 
solid line). The oscillation periods of the two frequency com- 
ponents are Ti = 1.5 and T2 = 2, thus the repetition time 
is Trep = [Ti,T2]lcm = 6. Noticc that the periodicity of the 
envelope (black dash line) is Tenv = 12, which is different from 
the repetition time Trep — 6. 



tributed^ polarization basis as described in Eq. (37) 



' COS c/ COS (p COS X — sm cp sm x ^ 
cos 6 sin (j) cos x + cos (j) sin x 
— sin cos X 

- cos 6 cos ^ sin X — sin cos x ^ 

- cos 6 sin (/) sin x + cos (/) cos x 
sin ^ sin x 



(B5) 



Appendix C: Repetitive Time 

A field composed of finite discrete frequencies, 

N 



E{t) = ^Ek cos (cjfet), 



(CI) 



/c=l 



repeats itself at the least common multiple (LCM) of all 
the periods of its frequency components, 



Trep = [TijT2j . . . ^ T^] LCM , 



(C2) 
(C3) 



27r 



where = — . An example of a two- frequency beat 
wave is given in Fig. [14| 



After the rotation, the uniformly distributed circle will span into 
a uniformly distributed spherical surface. 



Given the frequency spectrum of E{t)^ one can draw a 
relation between the repetition time Trep and the greatest 
common divider (GCD) of the frequencies. 



27T 



(a;i,cj2, • • • ,^n)gcd 



(04) 



The derivation of the relation in Eq. ( |C4| ) is the following. 
First, we factorize the sum of all the frequencies into two 
terms, 



+ 002 


+ ...4 


27r 


27r 








27T 



27T 



(C5) 



-(ni + n2 + . . . + un), 



[^1,^2, . . . ,Tn]lcm 
where Uk are positive integers, 

[^1,^2, . . . ,Tn]lcm 



rik 



Tk 



Now, it is true that 

(ni,n2, . . ■,nN)GCD = 1, 



(06) 
(07) 



because the reverse, (ni, n2 , . . . ^^1^)000 > 1, would lead 
to a contradiction to Eq. (C6). Therefore, one can con- 
clude that 



2lT 



(C8) 



[^1, ^2, • • • , Tn]lcm 



From Eq. ([C3]) and Eq. ([C8]), the relation in Eq. ( |C4 ) 

27r 



(09) 



is drawn. 

Since the simulation should only be carried through an 
integration time Tint < 'Trep to avoid repetitive solutions, 
in our case the choice of the integration time (Eq. (67)) 



27r 



(CIO) 



where Auj is the smallest frequency gap {Auj < Aujgcd)^ 
suffices our purpose. The frequency gap as a function of 
hz can be estimated using Eq. (33), 



Auj{k,) = uj{k,) — uj{k, — Atv) 
= c(3/t)i/3 



c(3«)i/3 ( 1 _ ^ 



1/3 



(Cll) 



Applying the sharp resonance condition (Eq. ([8|) to 
Eqs. (23) and (34), it can be further shown that Ati is 



much smaller than and hz c:^ hzo^ 
An _ f 3 \ A 

1^0 V^k -1/^0 



«1, 



(C12) 
(C13) 
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1 / (jJq \ 3 

where tzo = - ( — • Therefore, the size of Auj(hz) is ap- 
3 V c / 

proximately fixed within the sampled frequency range A, 
and the smallest frequency gap Au can be approximated 



as 



A.^'-^^^. (C14) 
3 f^o 
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